1 Background

Major questions:

  1. Did we induce CTA in red foxes?
  2. How long did it take to induce CTA?
  3. Did we induce CTA in ravens?

2 Data preparation

Firstly, we’re going to install and load our packages. We will use the pacman Package Management Tool which allows us to install and load subsequent packages in a condensed and efficient way.

# Install pacman package
install.packages('pacman')
# Install and load required packages
pacman::p_load(dplyr, ggmap, ggspatial, ggplot2, ggpubr, janitor, lme4, lmerTest, lwgeom, MuMIn, readxl, sf, sp, tidyr, tidyverse, viridis)

Secondly, we’ll read in our data and explore it.

# Read in data
data <- read_excel(path="data.xlsx", sheet="data") %>%
  clean_names()

Let’s view the data.

2.1 Bait take

Now, we’ll wrangle some of our data to get it ready for analyses and plotting.

data <- read_excel(path="data.xlsx", sheet="data") %>%
  clean_names() %>%
  mutate(date = as.Date(date)) %>%
  mutate(bait_take = as.numeric(bait_take)) %>%
  mutate(phase = as.factor(phase)) %>%
  drop_na(bait_take) %>%
  mutate(bait_yes = ifelse(bait_take == 1, "Yes", "No"))

Ellie, explain why you’re going to exclude Fridays here

# Create a column for day of the week from the date variable


# Remove Fridays from the dataframe

2.2 Weather

The weather data were acquired from the Australian Bureau of Meteorology weather station 70351.

weather <- read_excel(path="data.xlsx", sheet="weather") %>%
  clean_names() %>%
  mutate(date = as.Date(date))

# Combine bait take and weather data into a single df
data <- left_join(data, weather, by=c("date"))

2.3 Species

# Subset to actual instances of bait take by a vertebrate
data <- data %>%
  subset(species != "Lizard") %>%
  mutate(fox_visited = ifelse(fox_visit_1 != 0, "Yes", 
                       ifelse(fox_visit_2 !=0, "Yes",
                       ifelse(fox_visit_3 !=0, "Yes", 
                       ifelse(fox_visit_4 !=0, "Yes", "No"))))) %>%
  
  mutate(fox_visits = ifelse(fox_visit_4 != 0, 4, 
                      ifelse(fox_visit_3 !=0, 3,
                      ifelse(fox_visit_2 !=0, 2, 
                      ifelse(fox_visit_1 !=0, 1, 0))))) %>%

  mutate(raven_visited = ifelse(raven_visit_1 != 0, "Yes", 
                         ifelse(raven_visit_2 !=0, "Yes",
                         ifelse(raven_visit_3 !=0, "Yes", "No")))) %>%
  
  mutate(raven_visits = ifelse(raven_visit_3 !=0, 3, 
                        ifelse(raven_visit_2 !=0, 2,
                        ifelse(raven_visit_1 !=0, 1, 0))))

write.csv(data, "data processed.csv")

3 Analyses

3.1 Bait take

3.1.1 Plots

3.1.1.1 All species

Bait take by treatment

# Read in processed data
data <- read.csv("data processed.csv") %>%
  clean_names() %>%
  mutate(date = as.Date(date))

# Plot by treatment type
bait_treat <- ggplot(data, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("All species") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_treat)

Bait take by treatment and phase

# Read in processed data
data <- read.csv("data processed.csv") %>%
  clean_names() %>%
  mutate(date = as.Date(date))

# Plot by treatment type
bait_treat_phase <- ggplot(data, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("All species") +
  facet_wrap(~phase, scales="free") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_treat_phase)

Percentage bait taken by each species by phase

# Subset to only potential foxes
taken <- data %>%
  subset(bait_take != 0) %>%
  subset(species != "Not taken") %>%
  group_by(species, treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)

# Plot by treatment type
taken_species_phase_hist <- ggplot(data=taken, 
                                        aes(x = species, y = perc, 
                                            col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  facet_wrap(~phase) +
  ggtitle("All species") +
  xlab("Species") + 
  ylab("Percentage of taken baits (%)") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(taken_species_phase_hist)

Plot bait take by rainfall

# Plot by rainfall type
bait_rain <- ggplot(data, aes(x = rainfall_mm, 
                                          y = bait_take, 
                                          col=treatment, 
                                          fill=treatment)) +
  geom_smooth() + 
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("All species") +
  xlab("Rainfall") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")
  
# Display the plot
print(bait_rain)

3.1.1.2 Foxes

Percentage fox visits by treatment and phase

# Subset to only potential foxes
fox_stats <- data %>%
  subset(species == "Fox") %>%
  group_by(treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)
  
# Plot by treatment type
visits_fox_histogram <- ggplot(data=fox_stats, 
                              aes(y = perc, x = treatment, 
                                  col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  facet_wrap(~phase) +
  ggtitle("Potential fox visits (including unknowns)") +
  xlab("Species") + 
  ylab("Percentage of fox visits (%)") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(visits_fox_histogram)

# Subset to only potential foxes
foxes_only <- data %>%
  subset(species != "Raven") %>%
  subset(species != "Unknown")
  
# Plot by treatment type
bait_treat_phase <- ggplot(data=foxes_only, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("Foxes vs no take (excluding unknowns)") +
  facet_wrap(~phase, scales="free") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_treat_phase)

3.1.1.3 Ravens

Percentage fox visits by treatment and phase

# Subset to only potential foxes
raven_stats <- data %>%
  subset(species == "Raven") %>%
  group_by(treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)
  
# Plot by treatment type
visits_raven_hist <- ggplot(data=raven_stats, 
                              aes(y = perc, x = treatment, 
                                  col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  facet_wrap(~phase) +
  ggtitle("Potential raven visits (including unknowns)") +
  xlab("Species") + 
  ylab("Percentage of raven visits (%)") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(visits_raven_hist)

# Subset to only potential foxes
ravens_only <- data %>%
  subset(species != "Fox") %>%
  subset(species != "Unknown")
  
# Plot by treatment type
bait_raven_treat_phase <- ggplot(data=ravens_only, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("Ravens vs no take (excluding unknowns)") +
  facet_wrap(~phase, scales="free") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_raven_treat_phase)

Bait take by fox visits

# Plot bait take by species visits
bait_visits <- ggplot(data, mapping=aes(x = as.Date(date))) +
  geom_smooth(data, mapping=aes(y = as.numeric(bait_take)), 
              color = "red", fill = "red") +
  geom_smooth(data, mapping=aes(y = as.numeric(fox_visits)), 
              color = "darkorange", fill = "darkorange")  +
  geom_smooth(data, mapping=aes(y = as.numeric(raven_visits)), 
              color = "black", fill = "black")  +
  theme_minimal() + 
  facet_wrap(~treatment) +
  scale_y_continuous(name = "Bait take (%)", limits = c(0, 1), 
                     sec.axis = sec_axis(~., name = "Fox (orange) and raven (black) visits", 
                                         breaks = c(0, 1))) +
  xlab("Date")

# Display the plot
print(bait_visits)

3.1.2 Models

3.1.2.1 All species

Build generalised linear models for our research questions.

# Model for bait take by treatment
summary(glm(bait_take ~ treatment, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment, data = data)
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.66038    0.02223  29.705  < 2e-16 ***
## treatmentTreatment -0.10124    0.03164  -3.199  0.00142 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2357498)
## 
##     Null deviance: 224.02  on 941  degrees of freedom
## Residual deviance: 221.60  on 940  degrees of freedom
## AIC: 1316.1
## 
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and rainfall
summary(glm(bait_take ~ treatment + rainfall_mm, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment + rainfall_mm, data = data)
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.638871   0.023237  27.493  < 2e-16 ***
## treatmentTreatment -0.100640   0.031505  -3.194  0.00145 ** 
## rainfall_mm         0.006263   0.002060   3.040  0.00243 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2337005)
## 
##     Null deviance: 224.02  on 941  degrees of freedom
## Residual deviance: 219.44  on 939  degrees of freedom
## AIC: 1308.9
## 
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and phase
summary(glm(bait_take ~ treatment + phase, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment + phase, data = data)
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.69484    0.05679  12.234  < 2e-16 ***
## treatmentTreatment -0.10146    0.03165  -3.205  0.00139 ** 
## phase              -0.01787    0.02710  -0.659  0.50980    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2358916)
## 
##     Null deviance: 224.02  on 941  degrees of freedom
## Residual deviance: 221.50  on 939  degrees of freedom
## AIC: 1317.7
## 
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and phase
summary(glm(bait_take ~ treatment + site, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment + site, data = data)
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.901039   0.148702   6.059 2.01e-09 ***
## treatmentTreatment -0.102424   0.065102  -1.573 0.116011    
## siteF01             0.046330   0.173125   0.268 0.789062    
## siteF02             0.201385   0.170912   1.178 0.238991    
## siteF03            -0.014532   0.174168  -0.083 0.933523    
## siteF04            -0.157455   0.168687  -0.933 0.350858    
## siteF05            -0.134257   0.170471  -0.788 0.431161    
## siteF06            -0.051039   0.171986  -0.297 0.766718    
## siteF07             0.098961   0.170001   0.582 0.560632    
## siteF08            -0.003736   0.173070  -0.022 0.982781    
## siteF09             0.161306   0.167225   0.965 0.335005    
## siteF10            -0.025888   0.171773  -0.151 0.880239    
## siteF11            -0.003736   0.173070  -0.022 0.982781    
## siteF12            -0.032747   0.169541  -0.193 0.846886    
## siteF13            -0.466256   0.169131  -2.757 0.005956 ** 
## siteF14            -0.560520   0.172711  -3.245 0.001216 ** 
## siteF15            -0.585677   0.170318  -3.439 0.000611 ***
## siteF16            -0.374846   0.168687  -2.222 0.026524 *  
## siteF17            -0.310130   0.170001  -1.824 0.068444 .  
## siteF18            -0.480433   0.171773  -2.797 0.005270 ** 
## siteF19            -0.621452   0.171156  -3.631 0.000298 ***
## siteF20            -0.571342   0.171773  -3.326 0.000916 ***
## siteF21            -0.651039   0.171986  -3.785 0.000164 ***
## siteF22            -0.610447   0.170471  -3.581 0.000361 ***
## siteF23            -0.466256   0.169131  -2.757 0.005956 ** 
## siteF24            -0.578201   0.169541  -3.410 0.000678 ***
## siteF25            -0.453736   0.173070  -2.622 0.008898 ** 
## siteF26            -0.434979   0.171773  -2.532 0.011502 *  
## siteF27            -0.466256   0.169131  -2.757 0.005956 ** 
## siteF28             0.201385   0.170912   1.178 0.238991    
## siteF29             0.196729   0.171156   1.149 0.250693    
## siteF30             0.106147   0.172711   0.615 0.538981    
## siteF31            -0.232064   0.172070  -1.349 0.177788    
## siteF32             0.058162   0.169541   0.343 0.731635    
## siteF34            -0.084329   0.172711  -0.488 0.625480    
## siteF35             0.008052   0.170001   0.047 0.962232    
## siteF40            -0.666907   0.171156  -3.896 0.000105 ***
## siteF48            -0.751374   0.174168  -4.314 1.78e-05 ***
## siteF52             0.104719   0.167760   0.624 0.532646    
## siteF53             0.103414   0.168687   0.613 0.539996    
## siteF54            -0.537745   0.170912  -3.146 0.001708 ** 
## siteF55            -0.091515   0.170949  -0.535 0.592552    
## siteF56             0.008601   0.170471   0.050 0.959774    
## siteF57            -0.526493   0.167171  -3.149 0.001690 ** 
## siteF58            -0.810130   0.170001  -4.765 2.20e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1493372)
## 
##     Null deviance: 223.41  on 937  degrees of freedom
## Residual deviance: 133.36  on 893  degrees of freedom
##   (4 observations deleted due to missingness)
## AIC: 924.16
## 
## Number of Fisher Scoring iterations: 2

3.1.2.2 Foxes

# Do fox visits change with treatment or phase?
summary(glm(fox_visits ~ treatment + phase + 
              rainfall_mm + maximum_temperature_c, data=data))
## 
## Call:
## glm(formula = fox_visits ~ treatment + phase + rainfall_mm + 
##     maximum_temperature_c, data = data)
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           -0.181054   0.226700  -0.799   0.4247  
## treatmentTreatment     0.022190   0.050005   0.444   0.6573  
## phase                  0.036098   0.053793   0.671   0.5024  
## rainfall_mm            0.007035   0.003418   2.058   0.0398 *
## maximum_temperature_c  0.014803   0.009938   1.490   0.1367  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.5886676)
## 
##     Null deviance: 556.46  on 941  degrees of freedom
## Residual deviance: 551.58  on 937  degrees of freedom
## AIC: 2181.1
## 
## Number of Fisher Scoring iterations: 2

3.1.2.3 Ravens

# Model for visits against treatment and phase
summary(glm(raven_visits ~ treatment + phase + 
              rainfall_mm + maximum_temperature_c, data=data))
## 
## Call:
## glm(formula = raven_visits ~ treatment + phase + rainfall_mm + 
##     maximum_temperature_c, data = data)
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)            0.120731   0.101268   1.192    0.233
## treatmentTreatment     0.016744   0.022338   0.750    0.454
## phase                 -0.002641   0.024030  -0.110    0.913
## rainfall_mm           -0.002171   0.001527  -1.422    0.155
## maximum_temperature_c -0.001225   0.004439  -0.276    0.783
## 
## (Dispersion parameter for gaussian family taken to be 0.1174665)
## 
##     Null deviance: 110.37  on 941  degrees of freedom
## Residual deviance: 110.07  on 937  degrees of freedom
## AIC: 662.88
## 
## Number of Fisher Scoring iterations: 2

3.2 GIS

3.2.1 Map responses

Here we map our responses of bait_take, fox_visits, and raven_visits across Ginninderry Conservation Corridor.

# Use Google API to fetch a base map
# ggmap::register_google(key="[add API key here]", write=TRUE)
map <- get_map(location=c(lon=148.9811, lat=-35.2350),
               zoom=15, source="google", maptype="satellite", crop=FALSE)
# Read in shapefile
gis <- st_read("shapefiles/bait_stations_egg_cta.shp") %>%
  fortify() %>%
  clean_names() %>%
  rename(site=name)
## Reading layer `bait_stations_egg_cta' from data source 
##   `X:\Honours\R\shapefiles\bait_stations_egg_cta.shp' using driver `ESRI Shapefile'
## Simple feature collection with 43 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 148.9673 ymin: -35.24309 xmax: 148.9927 ymax: -35.22524
## Geodetic CRS:  GDA2020
# Map the bait stations
gin_map <- ggmap(map) + 
  geom_point(data=gis, size=2, aes(x=long, y=lat, col=treatment)) + 
  theme_minimal() +
  xlab("") + ylab("") + labs(col="Treatment")

# Display the plot
print(gin_map)

# Combine the bait take df with the GIS df
data <- left_join(data, gis, by="site")

# Generalised linear model
mod <- glm(bait_take ~ lat + long + elevation, data=data)
summary(mod)
## 
## Call:
## glm(formula = bait_take ~ lat + long + elevation, data = data)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.841e+03  7.170e+02  -6.753 2.56e-11 ***
## lat         -1.071e+01  3.202e+00  -3.345 0.000857 ***
## long         2.997e+01  4.978e+00   6.020 2.51e-09 ***
## elevation   -7.231e-05  6.421e-04  -0.113 0.910369    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1887321)
## 
##     Null deviance: 222.12  on 930  degrees of freedom
## Residual deviance: 174.95  on 927  degrees of freedom
##   (11 observations deleted due to missingness)
## AIC: 1095.7
## 
## Number of Fisher Scoring iterations: 2
# Map the bait stations
mod_map <- ggmap(map) + 
  geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
             aes(x=long, y=lat, fill=bait_take)) + 
  theme(axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.ticks.y=element_blank(), 
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(fill="#363842"), 
        strip.text=element_text(colour="white")) +
  facet_wrap(~phase) +
  scale_fill_viridis_c() +
  xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Bait take")

# Display the plot
print(mod_map)

# Map the fox visitation
mod_map <- ggmap(map) + 
  geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
             aes(x=long, y=lat, fill=fox_visits)) + 
  theme(axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.ticks.y=element_blank(), 
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(fill="#363842"), 
        strip.text=element_text(colour="white")) +
  facet_wrap(~phase) +
  scale_fill_viridis_c() +
  xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Fox visits")

# Display the plot
print(mod_map)

# Map the fox visitation
mod_map <- ggmap(map) + 
  geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
             aes(x=long, y=lat, fill=raven_visits)) + 
  theme(axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.ticks.y=element_blank(), 
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(fill="#363842"), 
        strip.text=element_text(colour="white")) +
  facet_wrap(~phase) +
  scale_fill_viridis_c() +
  xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Raven visits")

# Display the plot
print(mod_map)

3.2.2 Spatial features

Here we read in, project, and map the bait stations across Ginninderry Conservation Corridor (GCC) using the following shapefiles:

  1. Bait station points,
  2. GCC boundary,
  3. Roads,
  4. Construction area,
  5. Native vegetation areas,
  6. Pondage, and
  7. Watercourses.

What we need to do for each layer is check its inherent Coordinate Reference System (CRS, using st_transform() to find the last 4-digit number of its output), remind the layer of that number (using st_set_crs()), and then reproject it to the CRS we wish to work with (using st_transform("EPSG:4326")).

3.2.2.1 Calculate distance to features

  1. Bait station points
# Read in bait station points
stations <- st_read("shapefiles/bait_stations_egg_cta.shp") %>%
  as.data.frame() %>%
  fortify() %>%
  clean_names()
## Reading layer `bait_stations_egg_cta' from data source 
##   `X:\Honours\R\shapefiles\bait_stations_egg_cta.shp' using driver `ESRI Shapefile'
## Simple feature collection with 43 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 148.9673 ymin: -35.24309 xmax: 148.9927 ymax: -35.22524
## Geodetic CRS:  GDA2020
# Tranform to Spatial Points Dataframe
stations_spdf <- SpatialPointsDataFrame( 
  data.frame(stations$long, stations$lat), 
  stations, proj4string=CRS("EPSG:4326")) %>% 
  st_as_sf()
  1. GCC boundary
# Read in the layer
gcc_bou <- st_read("shapefiles/gcc_boundary_polygon.shp")
## Reading layer `gcc_boundary_polygon' from data source 
##   `X:\Honours\R\shapefiles\gcc_boundary_polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 148.9616 ymin: -35.24409 xmax: 148.9934 ymax: -35.21299
## Geodetic CRS:  WGS 84
# Remind the layer of its CRS then reproject to our desired CRS
gcc_bou <- gcc_bou %>%
  # Set projection using sf
  st_set_crs(st_crs(gcc_bou)) %>%
  st_transform("EPSG:4326")

# Calculate distance of bait stations to nearest edge of polygon
dist_bou <- st_geometry(obj = gcc_bou) %>%
  st_cast(to = 'MULTILINESTRING') %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_bou)
  1. Roads and tracks
# Read in the layer
gcc_roads <- st_read("shapefiles/i5516_roads.shp")
## Reading layer `i5516_roads' from data source 
##   `X:\Honours\R\shapefiles\i5516_roads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2775 features and 17 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 148.5 ymin: -36 xmax: 150 ymax: -35
## Geodetic CRS:  GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_roads <- gcc_roads %>%
  # Set projection using sf
  st_set_crs(st_crs(gcc_roads)) %>%
  st_transform("EPSG:4326")

# Convert linestrings to a single multilinestring
gcc_roads <- st_union(gcc_roads)

# Calculate distance of bait stations to nearest edge of polygon
dist_roads <- st_geometry(obj = gcc_roads) %>%
  st_cast(to = 'MULTILINESTRING') %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_roads)
  1. Construction area
# Read in the layer
construct <- st_read("shapefiles/gcc_construction_area.shp")
## Reading layer `gcc_construction_area' from data source 
##   `X:\Honours\R\shapefiles\gcc_construction_area.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 148.9725 ymin: -35.22814 xmax: 148.9903 ymax: -35.21997
## Geodetic CRS:  GDA2020
# Remind the layer of its CRS then reproject to our desired CRS
construct <- construct %>%
  # Change EPSG here to match the last number from st_crs() output
  st_set_crs(st_crs(construct)) %>%
  st_transform("EPSG:4326")

# Calculate distance of bait stations to nearest edge of polygon
dist_con <- st_geometry(obj = construct) %>%
  st_cast(to = 'MULTILINESTRING') %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_con)
  1. Native vegetation areas
# Read in the layer
gcc_veg <- st_read("shapefiles/i5516_nativevegetationareas.shp")
## Reading layer `i5516_nativevegetationareas' from data source 
##   `X:\Honours\R\shapefiles\i5516_nativevegetationareas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 313 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 148.5 ymin: -36 xmax: 150 ymax: -35
## Geodetic CRS:  GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_veg <- gcc_veg %>%
  # Change EPSG here to match the last number from st_crs() output
  st_set_crs(st_crs(gcc_veg)) %>%
  st_transform("EPSG:4326")

# Check for and fix invalid geometries
gcc_veg <- st_make_valid(gcc_veg)

# Convert multiple features to a single feature
gcc_veg <- st_union(gcc_veg)

# Calculate distance of bait stations to nearest edge of polygon
dist_veg <- st_geometry(obj = gcc_veg) %>%
  st_cast(to = 'MULTILINESTRING') %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_veg)
  1. Pondage
# Read in the layer
gcc_pond <- st_read("shapefiles/i5516_pondageareas.shp")
## Reading layer `i5516_pondageareas' from data source 
##   `X:\Honours\R\shapefiles\i5516_pondageareas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 149.0126 ymin: -35.31815 xmax: 149.5838 ymax: -35.05082
## Geodetic CRS:  GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_pond <- gcc_pond %>%
  # Change EPSG here to match the last number from st_crs() output
  st_set_crs(st_crs(gcc_pond)) %>%
  st_transform("EPSG:4326")

# Convert multiple features to a single feature
gcc_pond <- st_union(gcc_pond)

# Calculate distance of bait stations to nearest edge of polygon
dist_pond <- st_geometry(obj = gcc_pond) %>%
  st_cast(to = 'MULTILINESTRING') %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_pond)
  1. Watercourse lines
# Read in the layer
gcc_wcl <- st_read("shapefiles/i5516_watercourselines.shp")
## Reading layer `i5516_watercourselines' from data source 
##   `X:\Honours\R\shapefiles\i5516_watercourselines.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3627 features and 15 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 148.5 ymin: -36.00002 xmax: 150 ymax: -35
## Geodetic CRS:  GDA94
# Remind the layer of its CRS then reproject to our desired CRS
gcc_wcl <- gcc_wcl %>%
  # Change EPSG here to match the last number from st_crs() output
  st_set_crs("EPSG:4283") %>%
  st_transform("EPSG:4326")

# Convert linestrings to a single multilinestring
gcc_wcl <- st_union(gcc_wcl)

# Calculate distance of bait stations to nearest edge of polygon
dist_wcl <- st_geometry(obj = gcc_wcl) %>%
  st_cast(to = 'MULTILINESTRING') %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_wcl)
# Rename name to site
stations <- stations_spdf %>%
  rename(site = name)

# Join response df to spatial df
data <- left_join(data, stations, by="site") %>%
  filter(!is.na(bait_take))

3.2.2.2 Model relationship between distance to features

Here we model the relationship between our response variables (i.e., bait_take, fox_vists, and raven_visits) and the distance between the bait station and nearby spatial features (i.e., boundary [dist_bou], roads [dist_roads], construction site [dist_con], vegetation strip [dist_veg], pondage [dist_pond], watercourses [dist_wcl]).

# Subset to only variables in the subsequent models
data_mod <- data %>%
  select(bait_take, treatment, fox_visits, raven_visits, 
         dist_bou, dist_roads, dist_con, dist_veg, dist_pond, dist_wcl) %>%
  na.omit()

# Fit the model with glm
mod <- glm(bait_take ~ treatment + dist_bou + dist_roads + dist_con + 
              dist_veg + dist_pond + dist_wcl, 
            data = data_mod, family=binomial, na.action = "na.fail")

# Display summary statistics
summary(mod)
## 
## Call:
## glm(formula = bait_take ~ treatment + dist_bou + dist_roads + 
##     dist_con + dist_veg + dist_pond + dist_wcl, family = binomial, 
##     data = data_mod, na.action = "na.fail")
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -2.7788067  7.0046783  -0.397 0.691583    
## treatmentTreatment -0.9047092  0.1729551  -5.231 1.69e-07 ***
## dist_bou           -0.0033458  0.0008050  -4.156 3.23e-05 ***
## dist_roads          0.0030916  0.0005053   6.119 9.42e-10 ***
## dist_con            0.0032323  0.0009439   3.424 0.000616 ***
## dist_veg            0.0016249  0.0007593   2.140 0.032349 *  
## dist_pond          -0.0002584  0.0012654  -0.204 0.838209    
## dist_wcl            0.0003089  0.0006154   0.502 0.615655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1247.77  on 930  degrees of freedom
## Residual deviance:  966.61  on 923  degrees of freedom
## AIC: 982.61
## 
## Number of Fisher Scoring iterations: 4
# Since dist_pond and dist_wcl was non-significant, remove from model
mod <- glm(bait_take ~ treatment + dist_bou + dist_roads + dist_con + dist_veg, 
            data = data_mod, family=binomial, na.action = "na.fail")

# Display summary statistics
summary(mod)
## 
## Call:
## glm(formula = bait_take ~ treatment + dist_bou + dist_roads + 
##     dist_con + dist_veg, family = binomial, data = data_mod, 
##     na.action = "na.fail")
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -4.2406271  0.5364397  -7.905 2.68e-15 ***
## treatmentTreatment -0.8817769  0.1671349  -5.276 1.32e-07 ***
## dist_bou           -0.0031113  0.0006613  -4.704 2.55e-06 ***
## dist_roads          0.0030243  0.0004669   6.478 9.29e-11 ***
## dist_con            0.0031826  0.0003490   9.118  < 2e-16 ***
## dist_veg            0.0018967  0.0002137   8.877  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1247.8  on 930  degrees of freedom
## Residual deviance:  966.9  on 925  degrees of freedom
## AIC: 978.9
## 
## Number of Fisher Scoring iterations: 4
# Perform model selection using dredge() with AICc
dredge(mod, rank = "AICc") %>%
  filter(delta <2)
## Global model call: glm(formula = bait_take ~ treatment + dist_bou + dist_roads + 
##     dist_con + dist_veg, family = binomial, data = data_mod, 
##     na.action = "na.fail")
## ---
## Model selection table 
##       (Int)   dst_bou  dst_con  dst_rds  dst_veg trt df   logLik AICc delta
## [1,] -4.241 -0.003111 0.003183 0.003024 0.001897   +  6 -483.448  979     0
##      weight
## [1,]      1
## Models ranked by AICc(x)

Plot bait take against dist_bou, dist_con, dist_roads, dist_veg.

# Plot the distances against bait take
dist_bou_plot <- ggplot(data = data_mod) +
  geom_smooth(aes(x=dist_bou, y=bait_take), 
              col="purple", fill="purple") + 
  facet_wrap(~treatment) + 
  xlab("Distance to nearest boundary (m)") + 
  ylab("Bait take (proportion)") + 
  theme_minimal()

dist_con_plot <- ggplot(data = data_mod) +
  geom_smooth(aes(x=dist_con, y=bait_take), 
              col="orange", fill="orange") + 
  facet_wrap(~treatment) + 
  xlab("Distance to nearest part of construction site (m)") + 
  ylab("Bait take (proportion)") + 
  theme_minimal()

dist_road_plot <- ggplot(data = data_mod) +
  geom_smooth(aes(x=dist_roads, y=bait_take), 
              col="darkgrey", fill="darkgrey") + 
  facet_wrap(~treatment) + 
  xlab("Distance to nearest road (m)") + 
  ylab("Bait take (proportion)") + 
  theme_minimal()

dist_veg_plot <- ggplot(data = data_mod) +
  geom_smooth(aes(x=dist_veg, y=bait_take), 
              col="green", fill="green") +
  facet_wrap(~treatment) + 
  xlab("Distance to nearest vegetation patch (m)") + 
  ylab("Bait take (proportion)") + 
  theme_minimal()

dist_plots <- ggarrange(dist_bou_plot, dist_con_plot, 
                        dist_road_plot, dist_veg_plot, 
                        nrow=2, ncol=2)

# Display the plot
print(dist_plots)

Next, repeat these models and plots for fox visits and raven visits.

3.2.2.3 Map bait stations and features

Map these layers

# Map the site
gcc_map <- ggplot() + 
  geom_sf(data=gcc_roads, col="darkgrey") + 
  geom_sf(data=construct, fill="orange", alpha=0.3) + 
  geom_sf(data=gcc_veg, fill="green", alpha=0.3) + 
  geom_sf(data=gcc_pond, col="lightblue") + 
  geom_sf(data=gcc_wcl, col="blue") + 
  geom_sf(data=gcc_bou, fill="purple", alpha=0.3) + 
  geom_sf(data=stations_spdf, aes(col=treatment)) + 
  geom_text(data=stations_spdf, aes(x=long, y=lat, label=name), size=3, color="black") +
  # Limit map to GCC extent
  coord_sf(xlim = c(st_bbox(stations_spdf)[[1]]-0.005, 
                    st_bbox(stations_spdf)[[3]]+0.005), 
           ylim = c(st_bbox(stations_spdf)[[2]]-0.005, 
                    st_bbox(stations_spdf)[[4]]+0.005)) +
  theme_minimal() +
  scale_fill_manual(values=viridis(18, begin=0.9, end=0.1)) +
  xlab("") + ylab("") + labs(col="Treatment")

# Display the plot
print(gcc_map)

Session information